Análisis de los ingresos por remesas de México

Clase series de tiempo, primavera 2024

1) Definición

La primera actividad de cualquier análisis y pronostico es definir el problema

En el libro Forecasting: Principles and Practice, capítulo primero, describe que es un pronóstico, que se puede pronosticar y cómo definir el problema, y algunos casos de estudios.

Lee más en el capítulo https://otexts.com/fpp3/determining-what-to-forecast.html

Descripción del las remesas

Es el envío de dinero de aquellas personas que radican en otra nación a su país de origen. Por ejemplo, los envíos de dinero que realizan los mexicanos que radican en Estados Unidos y Canadá a sus familias que viven en México.

Con el surgimiento de nuevas herramientas tecnológicas es más sencillo realizar esta transacción y el número de usuarios que la hacen se ha incrementado en los últimos años.

También existen las remesas enviadas por personas que se dedican a los negocios y tienen que pagar sueldos a los empleados contratados en el exterior, por ejemplo. Otra transferencia muy común es para estudiantes que no pueden trabajar y necesitan mantenerse en otro país mientras estudian.

Las remesas son la segunda fuente de divisas más importante de México después de los ingresos por petróleo. Representan cerca de 3% del PIB, 50% de las exportaciones petroleras, 135% de la inversión extranjera directa y 189% de los ingresos de viajeros internacionales.

Los ingresos por remesas provenientes del exterior ascendieron a 4,760 millones de dólares en diciembre de 2021, lo que implicó un aumento anual de 30.4%.

Por su parte, las remesas enviadas por residentes en México al exterior mostraron un crecimiento anual de 1.8%, al alcanzar un nivel de 114 millones de dólares.

Con estos resultados, el superávit de la cuenta de remesas de México con el resto del mundo fue de 4,646 millones de dólares, mayor al de 3,538 millones de dólares que se presentó en diciembre de 2020.

Con series desestacionalizadas, en el último mes de 2021 los ingresos y egresos por remesas exhibieron retrocesos mensuales de 1.9 y 2.3%, en igual orden. Así, en diciembre de 2021 el superávit de la cuenta de remesas se situó en 4,600 millones de dólares, que se compara con el de 4,687 millones de dólares que se observó en noviembre de 2021.

Para la totalidad de 2021, el valor de los ingresos por remesas fue de 51,594 millones de dólares, monto superior al de 40,605 millones de dólares reportado en 2020 y que significó una expansión anual de 27.1%.

Durante 2021, el 98.9% del total de los ingresos por remesas se realizó a través de transferencias electrónicas, al ubicarse en 51,045 millones de dólares. Por su parte, las remesas efectuadas en efectivo y especie2 y las money orders representaron el 0.7 y 0.4% del monto total, respectivamente, al registrar niveles de 333 y 216 millones de dólares, en el mismo orden.

Los egresos por remesas sumaron 1,057 millones de dólares en 2021, cifra mayor a la de 899 millones de dólares observada en 2020 y que implicó un incremento anual de 17.6%.

Con estos resultados, el saldo superavitario de la cuenta de remesas en 2021 fue de 50,537 millones de dólares, cifra superior a la de 39,706 millones de dólares reportada en 2020 y que representó un avance anual de 27.3%.

Datos

Banco de México registra las remesas como parte de la balanza de pagos y tiene registro mensual desde 1995.

Los datos de Banxico pueden ser descargados en el siguiente enlace

Determinar que pronosticar

Se analizara el ingreso total de remesas, como se muestra en la captura siguiente.

Por que es importante analizar y pronosticar en el contexto social

El flujo de remesas internacionales tiene un efecto importante sobre la actividad económica en algunas entidades que se localizan en las regiones centrales y el sur del país. Esto condujo a que el incremento observado en las remesas haya inducido en 2016 un mayor crecimiento económico en estas regiones.

En efecto, como resultado de los patrones de gasto de los hogares, dichas entidades tienden a ser más propensas a concentrar el impacto de las remesas en las actividades comerciales y los servicios. Adicionalmente, los flujos de remesas contribuyen a suavizar el ingreso de los hogares ante choques sobre la actividad económica, en particular en el sector agropecuario.

Las remesas son la segunda fuente de divisas más importante de México. En los últimos años las remesas han alcanzado máximos históricos, además, han obtenido mayor atención por las recientes investigaciones que señalan cómo cárteles del narco lavan dinero vía remesas.

2) Datos

Importar

En las siguientes celdas se muestan los pasos para importar los datos desde mi espacio de trabajo. Por que es un archivo en formato Excel, la librería usada es readxl.

Lee más en el notebook introductorio visto en clase: tidyverse, sección Data science workflow, paso 1 Import Data.

También en el libro R for Data Science capítulo 20.

library(readxl)

#raw_data = read_xlsx(path = "/Users/danielnuno/GitHub/time_series_s2024/Analisis/remesas banxico.xlsx", sheet = "Hoja1", col_names = FALSE)
raw_data = read_xlsx(path = "~/GitHub/time_series_s2024/Analisis/remesas/remesas banxico.xlsx", sheet = "Hoja1", col_names = FALSE)

New names:

head(raw_data, 20)

Eliminar filas

Los datos no inician hasta la fila 18. Usando dplyr::slice se pueden eliminan las primeras 18 filas.

Lee más en el notebook introductorio visto en clase: tidyverse, sección Data science workflow, paso 3 understand your data.

También en el libro R for Data Science capítulo 27.

#raw_data = raw_data[-(1:18), ]
data = dplyr::slice(raw_data, -(1:18))
head(data, 2)

Renombra las columnas

Las dimensiones están nombradas como …1 y …2, por lo tanto es conveniente renombrarlas sin espacios ni caracteres especiales.

colnames(data) = c('date', 'value')
head(data, 2)

Cambia el tipo de dato de las columnas

Las columnas son chr pero date tiene que ser una fecha mensual y value son millones de dólares. Utiliza lubridate para hacer la transformación de la fecha o R básico.

Lee más en el notebook introductorio visto en clase: tidyverse, sección Data science workflow, paso 3 understand your data: data transformation.

También en el libro R for Data Science capítulo 17.

Transforma el tipo texto cómo numérico

data$date = as.numeric(data$date)

Transforma el número a fecha utilizando el argumento origin

El número representa un día en la serie que inicia el primero de enero de 1990. Por lo tanto utilizamos el argumento origin que indica el inicio de la serie.

data$date = as.Date(data$date, origin = "1899-12-30")

Crea una secuencía explicita la periodicidad mensual

data$date = seq(from = min(data$date), to = max(data$date), by = "1 month")
library(lubridate)
library(tsibble)
data$date = yearmonth(data$date)
head(data, 2)

Cambia el tipo de dato de la columna value

data$value = as.numeric(data$value)
head(data, 2)

Crea un tsibble

Usando tsibble::as_tsibble() asigna la fecha como index. En este caso key no es necesario pero decide si es importante para tu serie.

Lee más en el notebook introductorio visto en clase: tidyverse, sección Data science workflow, paso 2 Tidy Data. Y el ejemplo 1.

También la documentación tsibble y el libro Forecast: Principles and Practice capítulo 2.1

data = as_tsibble(data, index=date, regular=TRUE)
head(data, 2)

interval(data)

<interval[1]>
[1] 1M

Valores faltantes y factores

Identifica si tu serie de tiempo tiene valores faltantes o factores y has las transformaciones correspondientes con ayuda del libro R for data science - Transform.

3) Análisis

Visualización de la serie

Utiliza ggplot2 para gráficar la serie

Lee más en el notebook introductorio visto en clase: tidyverse, sección Data science workflow, paso 3 Understand your data - visualization.

Ve los ejemplos y explicaciones vistas en clase: ggplot2, sección Gráficas con ggplot2.

También el libro Forecast: Principles and Practice capítulo 2.2

Time plot

library(ggplot2)
library(feasts)
feasts::autoplot(data) + ggtitle('Remesas por mes') + ylab('Millones de dólares') + xlab('Fecha')

Plot variable not specified, automatically selected `.vars = value`

Patrones

Gráficamente presenta tendencia, estacionalidad y heteroscedasticidad. Parece presentar ciclos relacionados a la actividad económica.

Consulta el libro Forecast: Principles and Practice capítulo 2.3

Gráfica estacionales

Una gráfica estacional es lo mismo que un time plot pero graficada usando los periodos estacionales. En este caso la estacionalidad es mensual.

Consulta el libro Forecast: Principles and Practice capítulo 2.4

data %>% gg_season(value, labels = "both") +
    ggtitle('Remesas por año') + ylab('Millones de dólares') + xlab('Mes')

Utiliza plotly::ggplotly() hacer gráficos interactivos

suppressWarnings(library(plotly))

yearly_data_plot = data %>% gg_season(value, labels = "both") +
    ggtitle('Remesas por año') + ylab('Millones de dólares') + xlab('Mes')

ggplotly(yearly_data_plot)

Sub gráficas estacionales

Un gráfico alternativo que enfatiza los patrones estacionales es aquel en el que los datos de cada temporada se recopilan en minigráficos de tiempo separados.

Consulta el libro Forecast: Principles and Practice capítulo 2.5

subseries_plot = data %>% gg_subseries(value)
ggplotly(subseries_plot)

Gráfico de rezagos

Estos gráficos son gráficos de dispersión evaluados en los rezagos de la misma serie que miden la correlación entre periodos. En otras palabras, para esta serie de tiempo mensual, mide la relación de un mes en el tiempo contra un mes (n lags) del pasado.

Consulta el libro Forecast: Principles and Practice capítulo 2.7

lags_plots = data %>% filter(year(date) > 2018) %>% gg_lag(value, geom = "point", lags = 1:12) + labs(x = "lag(Remesa, k)")

suppressWarnings(ggplotly(lags_plots))

Autocorrelación

Así como la correlación mide el alcance de una relación lineal entre dos variables, la autocorrelación mide la relación lineal entre valores rezagados de una serie temporal.

Consulta el libro Forecast: Principles and Practice capítulo 2.8

data %>% ACF(value, lag_max = 12)

data %>% ACF(value, lag_max = 24) %>% autoplot() + labs(title='Remesas por año')

Cuando los datos tienen tendencia, las autocorrelaciones para desfases pequeños tienden a ser grandes y positivas porque las observaciones cercanas en el tiempo también lo son en valor. Por lo tanto, el ACF de una serie temporal con tendencia tiende a tener valores positivos que disminuyen lentamente a medida que aumentan los rezagos.

Cuando los datos son estacionales, las autocorrelaciones serán mayores para los rezagos estacionales (en múltiplos del período estacional) que para otros rezagos.

Sera importante realizar las mismas gráficas una vez que la serie sea estacional.

Además de lo descrito en la primera parte del documento, con las gráficas podemos observar lo siguiente:

Gráficamente presenta tendencia, estacionalidad. Parece presentar ciclos relacionados a la actividad económica.

es interesante observar con este gráfico que a partir del año 2004 los datos se ven separados, es decir, la amplitud de los datos se vuelve más grande en cada año después que pasa. Podemos concluir tres cosas con esto:

• Los primeros siete u ocho años de la serie hay crecimiento, pero a partir del 2003 el crecimiento es exponencial de cómo era antes.
• En esas fechas se volvió más fácil hacer transferencias electrónicas y por lo tanto contabilizarlo en la balanza de pagos.
• El crecimiento puede ser también por la solidaridad de los migrantes.

En el 2001, en respuesta a la “Sociedad para la prosperidad” entre México y Estados Unidos, el Banco de México y los Bancos de la Reserva Federal de los Estados Unidos acordaron estudiar la posibilidad de interconectar sus sistemas de pagos. La interconexión ha provisto a nuestras economías de un mecanismo eficiente para intercambiar pagos entre cuentas bancarias de ambos países.

En octubre de 2003, los Bancos de la Reserva Federal de los Estados Unidos y el Banco de México conectaron sus sistemas de pago para el envío de los pagos de los pensionados del Gobierno de los Estados Unidos que radican en México.

Posteriormente, desde el 2 de febrero de 2004, los usuarios de las instituciones financieras de los Estados Unidos suscritas a Directo a México pueden enviar pagos a cualquier cuenta bancaria en México.

Desde julio de 2005, este servicio de pagos se registró con el nombre Directo a México

Entonces para considerar como periodos que han afectado la serie de tiempo tenemos:

• 2003-10 facilidad de transferencias electrónicas
• 2007-12, 2009-06 recesión global por crisis hipotecaria
• 2020-02, 2020-04 recesión global por COVID-19
• Apartir del 2018 hay un cambio de tendencia y el crecimieto es acelerado.

Con los gráficos de rezagos podemos apreciar que muestra mayor correlación con el periodo inmediato anterior (1) y el mismo periodo del año pasado (12). Es decir, se espera que se envié una similar cantidad de dinero comparada a estos dos periodos del pasado (1 y 12), pero para el periodo 12 se espera que se envié más dinero.

Estadística descriptiva

Medidas de tendencia central

Calcula la media, moda y mediana. Estos estadisticos pueden ser de la totalidad de la serie o por año o por periodo estacional.

print(paste('fecha inicial', min(data$date)))

[1] "fecha inicial 1995 Jan"

print(paste('fecha final', max(data$date)))

[1] "fecha final 2023 Dec"

print(paste('observaciones', nrow(data)))

[1] "observaciones 348"

print(paste('existen', sum(is.na(data)), 'datos faltantes'))

[1] "existen 0 datos faltantes"

summary(data[, 'value'])

     value       
 Min.   : 248.1  
 1st Qu.: 797.8  
 Median :1933.9  
 Mean   :1993.3  
 3rd Qu.:2453.4  
 Max.   :5817.8  

Observamos que, por mes, el promedio de millones de dólares es de $1,993.2 millones de dólares, el mínimo fue de $248.1 millones de dólares, mientras que el máximo fue de $5,817.8 millones.

boxplot = data %>% 
            mutate(year = year(date)) %>% 
            ggplot(aes(x = as.factor(year), y = value)) + 
            geom_boxplot() + 
            xlab('Año') + 
            ylab('Remesas')

ggplotly(boxplot)

Medidas de dispersión

Estos estadísticos que nos dan información al respecto de la variabilidad o separación de los datos generalmente respecto a la media. Nos ayudan a comprender la distribución de los datos y evitan tomar conclusiones erróneas al comparar distintos grupos.

sd(data$value)

[1] 1313.988

var(data$value)

[1] 1726565

library(EnvStats)
kurtosis(data$value)

[1] 0.5389753

skewness(data$value)

[1] 0.8972091

shapiro.test(data$value)

    Shapiro-Wilk normality test

data:  data$value
W = 0.91084, p-value = 1.773e-13

library(ggExtra)
p <- ggplot(data, aes(x=date, y=value)) + 
        geom_hline(yintercept =1000) + 
        geom_hline(yintercept =3000) +
        geom_point() + 
        ggtitle('Remesas por mes') + ylab('Millones de dólares') + xlab('Fecha')

ggMarginal(p, type='histogram', margins = 'y')

histogram = ggplot(data, aes(x = value)) +
  geom_histogram( bins = 20, fill = "black", color = "black", alpha = 0.5) +
  labs(title = "Histograma",
       x = "Value",
       y = "Densidad")

ggplotly(histogram)

La distribución presenta un sesgo de simetría negativo a la izquierda y aunque los valores no son tan grandes, alguna transformación sería útil antes de hacer la descomposición de la serie y los pronósticos.

Valores atípicos

La detección de outliers es importante para afinar nuestro pronóstico y eliminar las observaciones atípicas. Existen diferentes formas de definir valores atípicos. En este caso, si los outliers son 1.5 veces el rango intercuartílico por encima y debajo del 1 y 3 cuartil, sí hay valores atípicos.

ttl_m_dlrs <- data %>% select('value')
ttl_m_dlrs <- as.numeric(unlist(ttl_m_dlrs[,1]))
summary(ttl_m_dlrs)[2] - 1.5*IQR(ttl_m_dlrs) >= summary(ttl_m_dlrs)[1]

1st Qu. 
  FALSE 

summary(ttl_m_dlrs)[5] + 1.5*IQR(ttl_m_dlrs) <= summary(ttl_m_dlrs)[6]

3rd Qu. 
   TRUE 

summary(ttl_m_dlrs)[2] - 3*IQR(ttl_m_dlrs) >= summary(ttl_m_dlrs)[1]

1st Qu. 
  FALSE 

summary(ttl_m_dlrs)[5] + 3*IQR(ttl_m_dlrs) <= summary(ttl_m_dlrs)[6]

3rd Qu. 
  FALSE 

library(tidyverse)

p <- data %>% as_tibble %>% group_by(years=year(date)) %>%
    summarise(remesas=sum(value)) %>%
    arrange(desc(years))%>%
    mutate(change = (remesas/lead(remesas) - 1) * 100) %>% 
    filter(years > 1995) %>% 
    filter(years < 2023)

mean_growth <- data %>% as_tibble %>% group_by(years=year(date)) %>%
                    summarise(remesas=sum(value)) %>%
                    arrange(desc(years))%>%
                    mutate(change = (remesas/lead(remesas) - 1) * 100) %>% 
                    filter(years > 1995) %>% 
                    filter(years < 2022) %>%
                    summarise(mean(change))

mean_growth <- mean_growth$`mean(change)`

ggplot(p, aes(x=years, y=change)) +
    geom_line() +
    geom_hline(yintercept=mean_growth) +
    geom_hline(yintercept=0) +
    ggtitle('Cambio porcentual por año') + ylab('%') + xlab('Mes')

El crecimiento promedio es de 11.48%.

4) Pronósticos base

Define los periodos de prueba y entrenamiento

Dependiendo de la cantidad de datos y la periodicidad podrías definir diferente cantidad de observaciones. La regla de pulgar es usar 80% de los datos para entrenar y el restante 20% para validar los resultados.

Para esta serie de tiempo 6 meses de validación son suficientes.

También en el libro Forecast: Principles and Practice capítulo 5

train <- data %>% select(value) %>% filter_index("1995 Jan" ~ "2023 Jun")
test <- data %>% select(value) %>% filter_index("2023 Jul" ~ "2023 Dec")
tstng_prds <- 6
frcst_prds <- 6

Seasonal Naive

En el libro Forecast: Principles and Practice capítulo 5, y notebook de clase pronósticos base se muestran diferentes métodos iniciales para hacer estimaciones. Estos métodos iniciales son ligados a la experiencia en la serie y el análisis exploratorio realizado con anterioridad.

Utilizando la funciónfable::SNAIVE para el primer pronóstico base debido a que hay una clara estacionalidad mensual.

library(fable)
models_fit <- train %>% 
    model(`Seasonal naive` = SNAIVE(value))
models_tst <- models_fit %>% forecast(h = tstng_prds)
snaive_plot <- models_tst %>% autoplot(filter_index(data, "2018 Jan" ~ .)) +
    ggtitle('Seasonal Naive') + ylab('Remesas') + xlab('Mes')

snaive_plot

Intervalos de predicción

models_tst

models_tst %>% hilo(level = c(80, 95))

Errores de pronóstico

accuracy(models_fit)

(models_fit %>% forecast(h = tstng_prds) %>% accuracy(test))

Diagnostico de resiuales

aug = augment(models_fit)
aug

aug %>% pull(.resid) %>% mean(na.rm = TRUE)

[1] 175.1076

aug %>% autoplot(.resid) + xlab("Mes") + ylab("") +
  ggtitle("Residuales del método seasonal naïve")

aug %>%
  ggplot(aes(x = .resid)) +
  geom_histogram() +
  ggtitle("Histograma de los residuales")

aug %>% ACF(.resid)

aug %>% ACF(.resid) %>% autoplot() + ggtitle("ACF of residuals")

train %>% 
  model(SNAIVE(value)) %>% 
  gg_tsresiduals()

Observamos que estas gráficas tienen un comportamiento muy distinto al Naïve:

• En la gráfica de los residuales vemos que se distingue claramente un patrón. De hecho, es el mismo patrón exactamente que siguen los datos originales, restándoles su media.

• La función de autocorrelación tiene un comportamiento típico de una caminata aleatoria. Por lo tanto, las autocorrelaciones son significativas.

• El histograma de los residuos muestra claramente que no se distribuyen de manera normal.

Test de Ljung-Box

Un test relacionado y que, generalmente, es más preciso es el test de Ljung-Box.

En este caso es igual: valores grandes de la prueba son indicios de que las autocorrelaciones no provienen de ruido blanco.

Entonces, la hipótesis nula de estas pruebas es que la serie en cuestión no está autocorrelacionada. En otras palabras, la H0 dice que la serie es ruido blanco. Si α es el nivel de significancia (el nivel máximo de error que estamos dispuestos a aceptar) y si el ¨p-value <α, entonces rechazamos H0, de lo contrario, no rechazamos la H0.

aug %>% features(.resid, ljung_box, lag=12, dof=0)

5) Descomposición

Componentes y descomposición STL

Consulta el libro Forecast: Principles and Practice capítulo 3.6

Para estudiar las características de tendencia y estacionalidad se descompone y graficar usando el método STL. La función per de periodicity es la más adecuada al parámetro de suavización estacional s.windows. Es el componente estacional es el más pequeño, siendo la tendencia más significativa.

Podemos observar en el componente estacional captura, que el inicio de año es el más flojito y que, a mediados del año es cuando se envían más remesas.

El competente de tendencia y remanentes observamos y confirmamos lo mencionado anteriormente. A partir del 2004 el envío de remesas, o el registro, aumenta considerablemente y el en 2008 y 2020 el remanente es muy significativo. También, a partir del 2020 el componente remanente se acentúa.

Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. J. (1990). STL: A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics, 6(1), 3–33. http://bit.ly/stl1990

stl_model = data %>% dplyr::select(value) %>% stl(s.window = 'per')
plot(stl_model,main = 'Descomposicón de la serie con STL')

Transformaciones y adjustes

Consulta el libro Forecast: Principles and Practice capítulo 3.1

Ajustar los datos históricos a menudo puede conducir a series temporales más simples. ejemplos de pueden ser ajustes de calendario, ajustes de población, ajustes de inflación y transformaciones matemáticas. El propósito de estos ajustes y transformaciones es simplificar los patrones en los datos históricos eliminando fuentes conocidas de variación o haciendo que el patrón sea más consistente en todo el conjunto de datos. Los patrones más simples suelen ser más fáciles de modelar y conducen a pronósticos más precisos.

Cómo ya se mencionó anteriormente con el histograma, y con observar las magnitudes de los valores al inicio y final de la serie es conveniente una transformación de potencia. El histograma y el sesgo también muestra sesgo negativo. Además, la prueba de normalidad SHAPIRO-WILK (alpha = 0.05)

qqnorm(data$value)
qqline(data$value)

Valor p es menor que Alpha entonces no tenemos una distribución normal. De hecho, estamos lejos de una distribución normal pero el sesgo no esta tan mal siendo 0.5. Queremos lo más acercado a 0 de sesgo.

Aunque se puede entender y usar la transformación box cox para transformar una distribución no normal a una normal, el principal objetivo es corregir sesgos y varianzas desiguales de la siguiente manera:

\[ y(\lambda) =\begin{cases} \frac{(y^\lambda - 1)}{\lambda}& \lambda \neq 0 ,\\ log(y)& \lambda = 0. \end{cases} \]

El lambda óptimo se encuentra cuando se encuentra la distribución con el máximo loglikelihood o con el coeficiente del gráfico de probabilidad de correlación o Shapiro-Wilk.

Bickel, P. J., & Doksum, K. A. (1981). An analysis of transformations revisited. Journal of the American Statistical Association, 76(374), 296–311. DOI

Box, G. E. P., & Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society. Series B, Statistical Methodology, 26(2), 211–252. DOI

bc <- EnvStats::boxcox(data$value, lambda=c(-2, 2), optimize=TRUE, objective.name='Log-Likelihood')
bc_data <- EnvStats::boxcoxTransform(data$value, bc$lambda)

skewness(bc_data)

[1] -0.07935633

shapiro.test(bc_data)

    Shapiro-Wilk normality test

data:  bc_data
W = 0.9521, p-value = 3.251e-09

qqnorm(bc_data)
qqline(bc_data)

data <- data %>% mutate('remesas_trn' = bc_data)

6) Pronósticos base con STL y tranformación matemática

STL Seasonal Naive

models_fit <- train %>% 
  model(stlf = decomposition_model(
    STL(value ~ trend(window = 12), robust = TRUE),
    NAIVE(season_adjust)
  ))
models_tst <- models_fit %>% forecast(h = tstng_prds)
mape_sn <- (models_fit %>% forecast(h = tstng_prds) %>% accuracy(test))$MAPE
snaive_plot <- models_tst %>% autoplot(filter_index(data, "2018 Jan" ~ .)) +
    ggtitle('STL') + ylab('Remesas') + xlab('Mes')

snaive_plot

?decomposition_model

models_fit <- train %>% 
  model(stlf = decomposition_model(
    STL(log(value) ~ trend(window = 12), robust = TRUE),
    NAIVE(season_adjust)
  ))
models_tst <- models_fit %>% forecast(h = tstng_prds)
mape_sn <- (models_fit %>% forecast(h = tstng_prds) %>% accuracy(test))$MAPE
snaive_plot <- models_tst %>% autoplot(filter_index(data, "2018 Jan" ~ .)) +
    ggtitle('Seasonal Naive') + ylab('Remesas') + xlab('Mes')

snaive_plot

models_fit <- train %>% 
  model(
    `Seasonal naive` = SNAIVE(value),
    stlf = decomposition_model(
    STL(value ~ trend(window = 12), robust = TRUE),
    NAIVE(season_adjust)),
    log_stlf = decomposition_model(
            STL(log(value) ~ trend(window = 12), robust = TRUE),
            NAIVE(season_adjust))
  )
models_tst <- models_fit %>% forecast(h = tstng_prds)
mape_sn <- (models_fit %>% forecast(h = tstng_prds) %>% accuracy(test))$MAPE
snaive_plot <- models_tst %>% autoplot(filter_index(data, "2018 Jan" ~ .), level = NULL) +
    ggtitle('Diferentes modelos') + ylab('Remesas') + xlab('Mes')

snaive_plot

models_tst

models_tst %>% hilo(level = c(80, 95))

accuracy(models_fit)

train %>% 
  model(decomposition_model(
            STL(log(value) ~ trend(window = 12), robust = TRUE),
            NAIVE(season_adjust))) %>% 
  gg_tsresiduals()

models_fit[3] %>% gg_tsresiduals()

7) Regresión Lineal

Datos

La hipótesis radica en la relación y crecimiento económico de Estados Unidos precede y esta correlacionado a el incremento de remesas. Por lo que el índice Brave-Butters-Kelley puede ser relevante para la estimación de las remesas.

The Brave-Butters-Kelley Indexes (BBKI) are a research project of the Federal Reserve Bank of Chicago. The BBK Coincident and Leading Indexes and Monthly GDP Growth for the U.S. are constructed from a collapsed dynamic factor analysis of a panel of 500 monthly measures of real economic activity and quarterly real GDP growth.

Los datos mensuales se pueden obtener de la FRED. https://fred.stlouisfed.org/series/BBKMGDP

Importar los datos usando tidyquant

library(tidyquant)
gdp_us = tq_get("BBKMGDP", get="economic.data", from = "1995-01-01")

Convertir la fecha y dataframe a tsibble

gdp_us$date = yearmonth(gdp_us$date)
gdp_us = gdp_us %>% select(date, price) %>% as_tsibble(index=date, regular = TRUE)
gdp_us

Convertir el cambio percentual a base 100 con la suma acumulada de los cambios porcentuales

gdp_us = gdp_us %>% mutate(gdp_base = cumsum(price))

gdp_us

Separa entre entrenamiento y prueba

train_gdp <- gdp_us %>% select(gdp_base) %>% filter_index("1995 Jan" ~ "2023 Jun")
test_gdp <- gdp_us %>% select(gdp_base) %>% filter_index("2023 Jul" ~ "2023 Dec")

Renombrar las columnas

train_gdp = add_column(train_gdp, train$value) #append remesas
colnames(train_gdp)[3] = "remesas" 
colnames(train_gdp)[1] = "gdp"
test_gdp = add_column(test_gdp, test$value)
colnames(test_gdp)[3] = "remesas"
colnames(test_gdp)[1] = "gdp"

Gráfica de ambas series a traves del tiempo

train_gdp |>
  pivot_longer(c(remesas, gdp), names_to="Series") |>
  autoplot(value) +
  labs(y = "value")

Gráfica de disperción (correlación) entre ambas series

train_gdp %>% ggplot(aes(x = gdp, y = log(remesas))) +
  labs(y = "Remesas",
       x = "GDP Base 1") +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

Ajuste de regresión y reporte del modelo

fit_lm <- train_gdp |>
  model(tslm = TSLM(log(remesas) ~ gdp))

report(fit_lm)

Series: remesas 
Model: TSLM 
Transformation: log(remesas) 

Residuals:
     Min       1Q   Median       3Q      Max 
-0.47518 -0.16670 -0.03232  0.17573  0.58386 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 5.736e+00  2.848e-02  201.44   <2e-16 ***
gdp         3.479e-03  5.678e-05   61.26   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.227 on 340 degrees of freedom
Multiple R-squared: 0.9169, Adjusted R-squared: 0.9167
F-statistic:  3753 on 1 and 340 DF, p-value: < 2.22e-16

Tabla aumentada del ajuste del modelo

augment(fit_lm)

Valores calculados

plot_lm = augment(fit_lm) |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = remesas, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))
ggplotly(plot_lm)

Zoom a la gráfica despues de 2018

augment(fit_lm) |> filter(year(date) > 2018) %>%
  ggplot(aes(x = date)) +
  geom_line(aes(y = remesas, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))

Analisis de residuales

Primera gráfica: errores en el tiempo para evaluar media y varianza. Segunda gráfica: autocorrelación de errores: errores en el pasado afectan el valor actual. tercera gráfica: histograma de errores para verificar la normalidad, media sea 0 y sesgo sea 0.

fit_lm %>% gg_tsresiduals()

augment(fit_lm) %>% features(.innov, ljung_box, lag=12)

Los residuales estan autocorrelacionados, no son homocedasticos y la media de los errores no es cero.

Gráfica de residuales contra valores ajustados.

También sirve para identificar outliers

augment(fit_lm) %>% ggplot(aes(x=.fitted, y=.resid)) + geom_point() + labs(x="Fitted", y="Residuals")

Los datos de ambas variables no son estacionarias, es decir tienen tendencía. Y puede ser que ningua tenga relación con la otra ni sea causal de la otra.

Regresiones con relaciones espurias puede que funcionen en el corto plazo pero dejen de funcionar en el futuro.

Pronósticos

test_gdp[c("gdp", "date")]

Gráfica de prónosticos

fc_lm = forecast(fit_lm, new_data = test_gdp[c("gdp", "date")])
data %>% autoplot(value) + autolayer(fc_lm)

Table de prónosticos

fc_lm

Otras variables predictoras

Tendencia

\[ y_{t} = \beta_{0} + \beta_{1}t + \epsilon_{t} \] dónde $ t = {1,2,3, …, n} $

TSLM utiliza el argumento trend() para indicar la tendencia.

Estacionalidad

Para estacionalidad diaria, calcula 6 variables:

\[ y_{t} = \beta_{0} + \beta_{1}lunes + \beta_{2}martes + \beta_{3}miercoles + \beta_{4}jueves + \beta_{5}viernes + \beta_{6}sabado + \epsilon_{t} \]

Para estacionalidad mensual, calcula 11 variables:

\[ y_{t} = \beta_{0} + \beta_{1}enero + \beta_{2}febrero + \beta_{3}marzo + \beta_{4}abril + \beta_{5}mayo + \beta_{6}junio + \beta_{7}julio + \beta_{8}agosto + \beta_{9}septiembre + \beta_{10}octubre + \beta_{11}noviembre + \epsilon_{t} \]

Para estacionalidad trimestral calcula 3 variables. Para estacionalidad diaria calcula 23 variables. Para estacionalidad semanal calcula 51 variables.

TSLM utiliza el argumento season() para indicar la estacionalidad.

Atipicos/Spike/Intervention variable

Con lo observado en las gráfica y lo comentado anteriormente, las fechas de importancia para esta serie de tiempo son la crisis hipotecaria y la crisis covid. Definir una columna de 0s, exepto entre las fechas indicadas (diciembre 2007 - junio 2009). Definir una columna de 0s, exepto entre las fechas indicadas (febrero 2020 - mayo 2020).

Crea una variable dummy de valores atípicos para tu serie de tiempo

train_gdp = train_gdp %>%
    mutate(crisish = if_else(date >= as.Date("2007-12-01") & 
                               date <= as.Date("2009-06-01"), 1, 0),
           crisisc = if_else(date >= as.Date("2020-02-01") & 
                               date <= as.Date("2020-05-01"), 1, 0)
           )

Warning: There were 4 warnings in `mutate()`.
The first warning was:
ℹ In argument: `crisish = if_else(...)`.
Caused by warning in `is_logical()`:
! Incompatible methods (">=.vctrs_vctr", ">=.Date") for ">="
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.

test_gdp = test_gdp %>%
    mutate(crisish = if_else(date >= as.Date("2007-12-01") & 
                               date <= as.Date("2009-06-01"), 1, 0),
           crisisc = if_else(date >= as.Date("2020-02-01") & 
                               date <= as.Date("2020-05-01"), 1, 0)
           )

Warning: There were 4 warnings in `mutate()`.
The first warning was:
ℹ In argument: `crisish = if_else(...)`.
Caused by warning in `is_logical()`:
! Incompatible methods (">=.vctrs_vctr", ">=.Date") for ">="
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.

Rezagos

Usa los valores pasados cómo predictor de X. Hipotesis: Las remesas (o el cambio porcentual) enviadas en el pasado influyen en el futuro. Revisa la autocorrelación.

\(x_{1}\) = remesas enviadas en el mes previo.

\(x_{2}\) = remesas enviadas hace dos meses.

\(x_{m}\) = remesas enviadas hace m meses.

El siguiente código muestra como calcular una variable con rezagos = 1

train_gdp$lag1 = c(NA, train_gdp$remesas[1:length(train_gdp$remesas)-1])
train_gdp

test_gdp$lag1 = c(train_gdp$remesas[length(train_gdp$remesas)], test_gdp$remesas[1:length(test_gdp$remesas)-1])
test_gdp

fourier

fit_lm <- train_gdp |>
  model(tslm = TSLM(log(remesas) ~ fourier(K=2) ))

report(fit_lm)

Series: remesas 
Model: TSLM 
Transformation: log(remesas) 

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7513 -0.6527  0.2549  0.4639  1.3582 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          7.31039    0.04258 171.692   <2e-16 ***
fourier(K = 2)C1_12 -0.09054    0.06021  -1.504    0.134    
fourier(K = 2)S1_12  0.02135    0.06022   0.355    0.723    
fourier(K = 2)C2_12 -0.03742    0.06021  -0.621    0.535    
fourier(K = 2)S2_12 -0.03407    0.06021  -0.566    0.572    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7873 on 337 degrees of freedom
Multiple R-squared: 0.009207,   Adjusted R-squared: -0.002553
F-statistic: 0.7829 on 4 and 337 DF, p-value: 0.53692

plot_lm = augment(fit_lm) |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = remesas, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))
ggplotly(plot_lm)

Fourier

fit_lm <- train_gdp |>
  model(tslm = TSLM(log(remesas) ~ trend() + fourier(K=6) ))

report(fit_lm)

Series: remesas 
Model: TSLM 
Transformation: log(remesas) 

Residuals:
     Min       1Q   Median       3Q      Max 
-0.51478 -0.21600 -0.07182  0.11848  0.65234 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          6.0477230  0.0313721 192.774  < 2e-16 ***
trend()              0.0073633  0.0001585  46.445  < 2e-16 ***
fourier(K = 6)C1_12 -0.0902607  0.0221349  -4.078 5.71e-05 ***
fourier(K = 6)S1_12  0.0212783  0.0221349   0.961   0.3371    
fourier(K = 6)C2_12 -0.0296631  0.0221354  -1.340   0.1811    
fourier(K = 6)S2_12 -0.0213924  0.0221366  -0.966   0.3346    
fourier(K = 6)C3_12 -0.0037863  0.0221349  -0.171   0.8643    
fourier(K = 6)S3_12 -0.0185609  0.0221349  -0.839   0.4023    
fourier(K = 6)C4_12 -0.0195143  0.0221354  -0.882   0.3786    
fourier(K = 6)S4_12 -0.0166023  0.0221351  -0.750   0.4538    
fourier(K = 6)C5_12 -0.0091343  0.0221349  -0.413   0.6801    
fourier(K = 6)S5_12 -0.0463210  0.0221349  -2.093   0.0371 *  
fourier(K = 6)C6_12 -0.0042256  0.0156519  -0.270   0.7873    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2894 on 329 degrees of freedom
Multiple R-squared: 0.8693, Adjusted R-squared: 0.8645
F-statistic: 182.4 on 12 and 329 DF, p-value: < 2.22e-16

plot_lm = augment(fit_lm) |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = remesas, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))
ggplotly(plot_lm)

Regresión lineal múltiple

fit_lm <- train_gdp |>
  model(tslm = TSLM(log(remesas) ~ trend() + season() + gdp + crisish + crisisc + fourier(K = 2)))

report(fit_lm)

Series: remesas 
Model: TSLM 
Transformation: log(remesas) 

Residuals:
     Min       1Q   Median       3Q      Max 
-0.54428 -0.13482  0.01035  0.15911  0.37922 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          5.4355637  0.0451541 120.378  < 2e-16 ***
trend()             -0.0042571  0.0006211  -6.854 3.61e-11 ***
season()year2        0.0170931  0.0496706   0.344 0.730971    
season()year3        0.1857560  0.0496679   3.740 0.000217 ***
season()year4        0.1730603  0.0496670   3.484 0.000561 ***
season()year5        0.2917365  0.0496664   5.874 1.05e-08 ***
season()year6        0.2336911  0.0495595   4.715 3.58e-06 ***
season()year7        0.2109345  0.0500205   4.217 3.21e-05 ***
season()year8        0.2199883  0.0500269   4.397 1.49e-05 ***
season()year9        0.1514828  0.0500320   3.028 0.002661 ** 
season()year10       0.1658213  0.0500363   3.314 0.001023 ** 
season()year11       0.0642423  0.0500380   1.284 0.200099    
season()year12       0.1289277  0.0500094   2.578 0.010374 *  
gdp                  0.0053772  0.0002823  19.046  < 2e-16 ***
crisish              0.2735281  0.0452035   6.051 3.95e-09 ***
crisisc              0.2360271  0.0993590   2.375 0.018103 *  
fourier(K = 2)C1_12         NA         NA      NA       NA    
fourier(K = 2)S1_12         NA         NA      NA       NA    
fourier(K = 2)C2_12         NA         NA      NA       NA    
fourier(K = 2)S2_12         NA         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1887 on 326 degrees of freedom
Multiple R-squared: 0.9449, Adjusted R-squared: 0.9424
F-statistic:   373 on 15 and 326 DF, p-value: < 2.22e-16

Selección de variables

glance(fit_lm) |>
  select(adj_r_squared, CV, AIC, AICc, BIC)

fit_lm <- train_gdp |>
  model(tslm = TSLM(log(remesas) ~ trend() + season() + gdp + crisish + crisisc + lag1))

report(fit_lm)

Series: remesas 
Model: TSLM 
Transformation: log(remesas) 

Residuals:
     Min       1Q   Median       3Q      Max 
-0.48499 -0.11092 -0.01781  0.15005  0.36325 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     5.472e+00  4.303e-02 127.178  < 2e-16 ***
trend()        -5.233e-03  5.851e-04  -8.943  < 2e-16 ***
season()year2   6.024e-02  4.625e-02   1.303 0.193666    
season()year3   2.212e-01  4.615e-02   4.793 2.50e-06 ***
season()year4   1.550e-01  4.607e-02   3.364 0.000862 ***
season()year5   2.879e-01  4.599e-02   6.259 1.23e-09 ***
season()year6   1.916e-01  4.625e-02   4.142 4.40e-05 ***
season()year7   1.923e-01  4.639e-02   4.146 4.33e-05 ***
season()year8   2.078e-01  4.634e-02   4.484 1.02e-05 ***
season()year9   1.347e-01  4.637e-02   2.905 0.003921 ** 
season()year10  1.672e-01  4.630e-02   3.611 0.000353 ***
season()year11  5.491e-02  4.633e-02   1.185 0.236823    
season()year12  1.499e-01  4.632e-02   3.236 0.001338 ** 
gdp             5.007e-03  2.656e-04  18.853  < 2e-16 ***
crisish         2.375e-01  4.172e-02   5.692 2.81e-08 ***
crisisc         2.303e-01  9.117e-02   2.526 0.012006 *  
lag1            1.573e-04  1.984e-05   7.929 3.60e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1731 on 324 degrees of freedom
Multiple R-squared: 0.9532, Adjusted R-squared: 0.9509
F-statistic: 412.8 on 16 and 324 DF, p-value: < 2.22e-16

glance(fit_lm) |>
  select(adj_r_squared, CV, AIC, AICc, BIC)

Análisis de residuales

fit_lm %>% gg_tsresiduals()

8) Suavización exponencial

Simple

fit_es <- train |>
  model(ETS(value ~ error("A") + trend("N") + season("N")))

report(fit_es)

Series: value 
Model: ETS(A,N,N) 
  Smoothing parameters:
    alpha = 0.5842231 

  Initial states:
     l[0]
 261.8688

  sigma^2:  50869.13

     AIC     AICc      BIC 
5705.757 5705.828 5717.262 

Alpha es de 0.58, lo qué indica que la última observación no es la más importante y las pasadas observaciones tienen un peso significativo

AICc es 5705

plot_lm = augment(fit_es) |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))
ggplotly(plot_lm)

Con tendencia

fit_es <- train |>
  model(AAN = ETS(value ~ error("A") + trend("A") + season("N")))

report(fit_es)

Series: value 
Model: ETS(A,A,N) 
  Smoothing parameters:
    alpha = 0.5587472 
    beta  = 0.0001000121 

  Initial states:
     l[0]     b[0]
 270.2258 15.16794

  sigma^2:  50455.56

     AIC     AICc      BIC 
5704.948 5705.126 5724.122 

plot_lm = augment(fit_es) |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))
ggplotly(plot_lm)

fit_es <- train |>
  model(Damped = ETS(value ~ error("A") + trend("Ad") + season("N")))

report(fit_es)

Series: value 
Model: ETS(A,Ad,N) 
  Smoothing parameters:
    alpha = 0.5691449 
    beta  = 0.004834991 
    phi   = 0.9799996 

  Initial states:
     l[0]     b[0]
 269.8862 9.707265

  sigma^2:  51170.91

     AIC     AICc      BIC 
5710.749 5711.000 5733.758 

plot_lm = augment(fit_es) |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))
ggplotly(plot_lm)

Con estacionalidad

fit_es <- train |>
  model(AAA = ETS(value ~ error("A") + trend("A") + season("A")))

report(fit_es)

Series: value 
Model: ETS(A,A,A) 
  Smoothing parameters:
    alpha = 0.2878418 
    beta  = 0.01257488 
    gamma = 0.3483385 

  Initial states:
     l[0]     b[0]      s[0]     s[-1]   s[-2]     s[-3]    s[-4]    s[-5]    s[-6]    s[-7]     s[-8]    s[-9]
 391.0324 2.686417 -5.107151 -123.8252 71.0204 -5.325298 107.7276 90.86807 148.2313 194.8573 -67.60189 59.09909
    s[-10]    s[-11]
 -244.5589 -225.3853

  sigma^2:  22715.69

     AIC     AICc      BIC 
5443.656 5445.545 5508.848 

plot_lm = augment(fit_es) |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))
ggplotly(plot_lm)

train

fit_es <- train |>
  model(AAA = ETS(value ~ error("M") + trend("A") + season("M")))

report(fit_es)

Series: value 
Model: ETS(M,A,M) 
  Smoothing parameters:
    alpha = 0.4929797 
    beta  = 0.01305041 
    gamma = 0.1295237 

  Initial states:
     l[0]     b[0]      s[0]     s[-1]     s[-2]    s[-3]    s[-4]    s[-5]   s[-6]    s[-7]    s[-8]     s[-9]
 291.2129 3.110346 0.9814982 0.8821721 0.9925036 1.001262 1.081113 1.054762 1.07494 1.147514 1.014973 0.9830665
    s[-10]    s[-11]
 0.8696745 0.9165203

  sigma^2:  0.0048

     AIC     AICc      BIC 
5186.587 5188.476 5251.779 

plot_lm = augment(fit_es) |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = "reales")) +
  geom_line(aes(y = .fitted, colour = "ajustados")) +
  labs(y = NULL,
    title = "remesas"
  ) +
  guides(colour = guide_legend(title = NULL))
ggplotly(plot_lm)

Selección del modelo

fit_es <- train |>
  model(ETS(value))

report(fit_es)

Series: value 
Model: ETS(M,A,M) 
  Smoothing parameters:
    alpha = 0.4929797 
    beta  = 0.01305041 
    gamma = 0.1295237 

  Initial states:
     l[0]     b[0]      s[0]     s[-1]     s[-2]    s[-3]    s[-4]    s[-5]   s[-6]    s[-7]    s[-8]     s[-9]
 291.2129 3.110346 0.9814982 0.8821721 0.9925036 1.001262 1.081113 1.054762 1.07494 1.147514 1.014973 0.9830665
    s[-10]    s[-11]
 0.8696745 0.9165203

  sigma^2:  0.0048

     AIC     AICc      BIC 
5186.587 5188.476 5251.779 

components(fit_es) |>
  autoplot()

Pronostico, intervalo de predicción y análisis de residuales

frcst = fit_es %>% forecast(h = tstng_prds)

fit_es %>%
  forecast(h = tstng_prds) %>%
  autoplot(train)

frcst

frcst %>% hilo(level = c(80, 95))

accuracy(fit_es)

accuracy(fit_lm)

fit_es %>% gg_tsresiduals()